home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_2 / a2_6.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.8 KB  |  179 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 2.6 (Secant Method).
  9. % Section    2.4,    Newton-Raphson and Secant Methods, Page 85
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program implements the secant method.
  13.  
  14. %    Define and store the function f(x) in the M-file  f.m  
  15. % function y = f(x)
  16. % y = x.^3 - 3.*x + 2;
  17.  
  18. delete f.m
  19. diary f.m; disp('function y = f(x)');...
  20.            disp('y = x.^3 - 3.*x + 2;');...
  21. diary off;
  22.  
  23. % Remark. f.m and secant.m are used for Algorithm 2.6
  24. f(0); % Test for file f.m
  25. pause    % Press any key to see the graph y = f(x).
  26.  
  27. clc; clg
  28. a = -2.5;
  29. b = 2.5;
  30. c = -5;
  31. d = 5;
  32. h = (b-a)/150;
  33. X = a:h:b;
  34. Y = f(X);
  35. axis([a b c d]);...
  36. plot(X,Y,'-g');...
  37. hold on;...
  38. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  39. xlabel('x');...
  40. ylabel('y');...
  41. title('Graph of y = f(x).');...
  42. grid;...
  43. axis;...
  44. hold off;...
  45. shg; pause    % Press any key to perform secant iteration.
  46.  
  47. clc;
  48. %    Place the starting value in  p0  and  p1
  49.  
  50. %    Place the abscissa tolerance in  delta
  51.  
  52. %    Place the ordinate tolerance in  epsilon
  53.  
  54. %    Place the maximum number of iterations in  max1
  55.  
  56. p0 = -2.6;
  57. p1 = -2.4;
  58. delta = 1e-12;
  59. epsilon = 1e-12;
  60. max1  = 50;
  61.  
  62. [p1,y1,err,P] = secant('f',p0,p1,delta,epsilon,max1);
  63.  
  64. pause    % Press any key for the list of iterations.
  65.  
  66. clc; clg;
  67. a = -2.7;
  68. b = -1.9;
  69. c = -10;
  70. d = 2;
  71. h = (b-a)/150;
  72. X = a:h:b;
  73. Y = f(X);
  74. max1 = length(P);
  75. clear Vx Vy
  76. for i = 2:max1,
  77.   k1 = 3*i-2;
  78.   k2 = 3*i-1;
  79.   k3 = 3*i;
  80.   Vx(k1) = P(i);
  81.   Vy(k1) = 0;
  82.   Vx(k2) = P(i);
  83.   Vy(k2) = f(P(i));
  84.   Vx(k3) = P(i-1);
  85.   Vy(k3) = f(P(i-1));
  86. end
  87. Z = zeros(1,length(P));
  88. axis([a b c d]);...
  89. plot(X,Y,'-g',Vx,Vy,'-r',P,Z,'or');...
  90. hold on;...
  91. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  92. xlabel('x');...
  93. ylabel('y');...
  94. title('Graphical analysis for the secant method.');...
  95. grid;...
  96. axis;...
  97. hold off;...
  98. shg; pause    % Press any key to continue.
  99.  
  100. J = 1:(max1);
  101. Yp = f(P);
  102. points = [J;P;Yp];
  103. Mx1 = 'Iterations for the secant method.';
  104. Mx2 = '     k                  p(k)               f(p(k))';
  105. Mx3 = 'The solution is:';
  106. Mx4 = 'The error estimate for p is  ± ';
  107. clc,echo off, diary output,...
  108. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  109. disp('Iteration converged with order 1.618 to the simple root.'),...
  110. disp(''),disp(Mx3),disp(''),disp('p = '),...
  111. disp(p1),disp(''),disp('f(p) = '),disp(y1),...
  112. disp([Mx4,num2str(err)]),diary off,echo on
  113.  
  114. pause    % Press any key to perform secant iteration.
  115.  
  116. clc;
  117. %    Place the starting value in  p0  and  p1
  118.  
  119. %    Place the abscissa tolerance in  delta
  120.  
  121. %    Place the ordinate tolerance in  epsilon
  122.  
  123. %    Place the maximum number of iterations in  max1
  124.  
  125. p0 = 1.6;
  126. p1 = 1.4;
  127. delta = 1e-12;
  128. epsilon = 1e-12;
  129. max1  = 30;
  130.  
  131. [p1,y1,err,P] = secant('f',p0,p1,delta,epsilon,max1);
  132.  
  133. pause    % Press any key for the list of iterations.
  134.  
  135. clg;
  136. a = 0.95;
  137. b = 1.6;
  138. c = -0.1;
  139. d = 1.2;
  140. h = (b-a)/150;
  141. X = a:h:b;
  142. Y = f(X);
  143. max1 = length(P);
  144. clear Vx Vy
  145. for i = 2:max1,
  146.   k1 = 3*i-2;
  147.   k2 = 3*i-1;
  148.   k3 = 3*i;
  149.   Vx(k1) = P(i);
  150.   Vy(k1) = 0;
  151.   Vx(k2) = P(i);
  152.   Vy(k2) = f(P(i));
  153.   Vx(k3) = P(i-1);
  154.   Vy(k3) = f(P(i-1));
  155. end
  156. Z = zeros(1,length(P));
  157. axis([a b c d]);...
  158. plot(X,Y,'-g',Vx,Vy,'-r',P,Z,'or');...
  159. hold on;...
  160. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  161. xlabel('x');...
  162. ylabel('y');...
  163. title('Graphical analysis for the secant method.');...
  164. grid;...
  165. axis;...
  166. hold off;...
  167. shg; pause    % Press any key to continue.
  168.  
  169. J = 1:(max1);
  170. Yp = f(P);  % User defined function.
  171. points = [J;P;Yp];
  172. clc,echo off, diary output,...
  173. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  174. disp('Iteration is converging linearly to the double root.'),...
  175. disp(''),disp(Mx3),disp(''),disp('p = '),...
  176. disp(p1),disp(''),disp('f(p) = '),disp(y1),...
  177. disp([Mx4,num2str(err)]),diary off,echo on
  178.  
  179.